Chapter 9 Mixture of Gaussians and EM

如果我们定义了一个观测变量和隐变量的联合分布，那么观测变量的分布就可以通过对隐变量边缘化得到，也就是把隐变量积分掉。因此，可以把难以求解的观测变量的边缘分布扩展至一种相对容易处理的方式，即扩展到由观测变量和隐变量构成的联合分布，潜变量的引入可以把复杂的分布拆解成简单的分布。

[TOC]

混合模型一方面可以用于表达更复杂的分布，另一方面可以用于聚类。我们从聚类问题引入混合模型，首先讲经典的非概率聚类算法$[Math Processing Error]K$-means 算法，在这个算法中，隐变量可以解释为一个样本点对所属类别的分配系数。对于带有隐变量的模型，一种经典的最大似然估计器是expectation-maximization(EM)算法，我们首先用混合高斯分布来引入EM算法，然后用严谨的方式推广到一般型隐变量视角。可以发现，$[Math Processing Error]K$-means其实EM算法应用在高斯混合模型(GMM)上的一种非概率简化。

GMM广泛应用于数据挖掘、模式识别、机器学习和统计分析等领域，在大部分应用场景下，GMM的参数都是由最大似然得到的，尤其是EM法。当然，最大似然法是由许多显著缺陷的，后面章节会用变分推断来解决这个问题。相当于对EM算法做了一些扩充，变分推断解决了MLE的一些本质缺陷，并且可以自动设置混合模型中组件的数量。

9.1 $K$-means Clustering

首先考虑多维空间的聚类问题，假设数据集$\left\{\mathbf{x}_{1}, \ldots, \mathbf{x}_{N}\right\}$由$N$个$D$维欧式空间上的变量$\mathbf{x}$构成，我们的目标是把这些数据划分为$K$个聚集，首先假定$K$的数值是已知的，直觉上，分到一个簇的点，应该对这个簇中心的距离，小于对其他簇中心的距离，定义$D$维向量$\pmb{\mu}_{k}$表示第$k$个簇的中心坐标，其中$k=1,…,K$，那么我们聚类的目的就是找到每个点所属的簇，使得所有点到其簇中心的平方距离之和最小。

为了方便描述一个点对于一个簇的所属情况，我们引入一个二进制变量$r_{n k} \in\{0,1\}$，表示第$n$个点对第$k$个簇的归属情况，如果归属，则值为$r_{n k} =1$，如果不归属，则$r_{n k} =0$。这种编码方式称作1-of-$K$编码。此时，聚类问题的目标函数可以定义如下，也叫distortion measure:

$$\begin{equation} J=\sum_{n=1}^{N} \sum_{k=1}^{K} r_{n k}\left\|\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right\|^{2} \end{equation}$$

聚类目标为求解使得$J$达到最小时的$r_{n k}$和$\pmb{\mu}_{k}$，这个问题可以通过迭代优化求解，首先初始化$\pmb{\mu}_{k}$并固定，然后最小化$J$优化$r_{n k}$，然后再固定$r_{n k}$优化$\pmb{\mu}_{k}$，如此迭代直到收敛。在EM算法框架下，优化$r_{n k}$的过程就叫做Step Expectation，优化$\pmb{\mu}_{k}$的过程就叫做Step Maximization。当然，解决$K$-means问题本身不需要EM这么复杂的方式，这里只是抛砖引玉。

如果有了聚类中心$\pmb{\mu}_{k}$，每个样本的归属情况是可以直接计算出来的，样本最近的簇中心即为归属簇，用公式表示为：

$$r_{n k}=\left\{\begin{array}{ll}{1} & {\text { if } k=\arg \min _{j}\left\|\mathbf{x}_{n}-\boldsymbol{\mu}_{j}\right\|^{2}} \\ {0} & {\text { otherwise }}\end{array}\right.$$

得到归属情况$r_{n k}$后，就可以通过最小化$J$来求解下一轮的聚类中心$\pmb{\mu}_{k}$，通过对$J$求导，再使导数为零，可以得到聚类中心即为簇中所有点的平均值：

$$\boldsymbol{\mu}_{k}=\frac{\sum_{n} r_{n k} \mathbf{x}_{n}}{\sum_{n} r_{n k}}$$

以上即为$K$-means方法的两个迭代步骤，直到点的归属情况不变，就算收敛了，此算法的收敛过程是经过证明的，但是也可能收敛到局部最优化。事实上，$K$-means方法可以用于GMM的参数初始化。

注意，上面$\pmb{\mu}_{k}$的计算也可以用梯度下降法。距离函数也可不用平方误差，可以换成任意距离函数，如下：

$$\widetilde{J}=\sum_{n=1}^{N} \sum_{k=1}^{K} r_{n k} \mathcal{V}\left(\mathbf{x}_{n}, \boldsymbol{\mu}_{k}\right)$$

$K$-means算法是一种硬聚类方法，接下来引入一种基于概率的软聚类方法。

9.2 Mixture of Gaussians

从分布角度来看，混合高斯模型就是多个高斯分布的线性叠加，可以用于描述更加复杂的分布情况，并没有什么玄乎的。我们现在来考虑离散隐变量下的混合高斯模型，从而引入到EM算法。

首先，混合高斯模型的概率密度函数可以表示为：

$$p(\mathbf{x})=\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\mathbf{x} | \boldsymbol{\mu}_{k}, \mathbf{\Sigma}_{k}\right)$$

对于变量$\mathbf{x}$，我们引入一个1-of-$K$表示的归属变量 $\mathbf{z}$ ，是个向量，只有一个元素为1，其余都为零，也就是满足$\sum_{k} z_{k}=1$，因此向量 $\mathbf{z}$ 共有 $K$ 种状态。此时我们可以定义联合分布$p(\mathbf{x}, \mathbf{z})$，边缘分布 $p( \mathbf{z})$ 和条件分布 $p(\mathbf{x}| \mathbf{z})$， 边缘分布的概率，可以由GMM的混合系数$\pi_{k}$来表示：

$$p\left(z_{k}=1\right)=\pi_{k} ， p\left(z_{k}=0\right)=1-\pi_{k}$$

注意，混合参数满足$0 \leqslant \pi_{k} \leqslant 1$ 且 $\sum_{k=1}^{K} \pi_{k}=1$。由于向量中只有一个量为1，其他量都是0，那么 $p(\mathbf{z})$ 表示为：

$$p(\mathbf{z})=\prod_{k=1}^{K} \pi_{k}^{z_{k}}$$

式子中的指数项目只有一个位置是1，其他都是0。事件 $\mathbf{z}$ 发生的意思就是 $z_{k}=1$ ，此时样本只属于某一个高斯分布，于是条件概率即为：

$$p(\mathbf{x} | \mathbf{z})=\prod_{k=1}^{K} \mathcal{N}\left(\mathbf{x} | \boldsymbol{\mu}_{k}, \mathbf{\Sigma}_{k}\right)^{z_{k}}$$

可以看出，指数项只有一个是1，其他都是0，因此就是一个单元高斯分布而已，联合分布表示为边缘分布 $p( \mathbf{z})p(\mathbf{x}| \mathbf{z})$，然后再把 $\mathbf{z}$ 边缘化(求和或者积分)，就得到了$p( \mathbf{x})$，如下，向量 $\mathbf{z}$ 共有 $K$ 种状态，因此变化化也就是对 $K$ 种状态下的联合分布求和，注意这里没有指数项：

$$p(\mathbf{x})=\sum_{\mathbf{z}} p(\mathbf{x,z}) =\sum_{\mathbf{z}} p(\mathbf{z})p(\mathbf{x} | \mathbf{z})=\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\mathbf{x} | \boldsymbol{\mu}_{k}, \mathbf{\Sigma}_{k}\right)$$

请注意，对于每个观测样本 $\mathbf{x}_{n}$ 都有一个对应的隐变量 $\mathbf{z}_{n}$ 。

可以看到，式子(10)和式子(5)是相同的，也就是说，我们现在已经把混合高斯模型转换成了包含隐变量的形式，到目前位置这么做好像没有任何意义，往下看。

现在我们把求解边缘分布 $p(\mathbf{x})$ 的问题转化成了求解联合分布 $p(\mathbf{x,z})$的问题，这在后面会显著简化问题。

这里再引入一个非常重要的概率，即给定观测值下属于某个高斯分布的概率，也就是后验概率，很自然，后验概率 $\gamma(z_{k})$ 即为观测值下为所有分布在此观测位置的累计加权的比例值，这个概率是个复合因素，既有归属系数，又有高斯分布的概率：

$$\begin{aligned} \gamma\left(z_{k}\right) \equiv p\left(z_{k}=1 | \mathbf{x}\right)=& \frac{p\left(z_{k}=1\right) p\left(\mathbf{x} | z_{k}=1\right)}{\sum_{j=1}^{K} p\left(z_{j}=1\right) p\left(\mathbf{x} | z_{j}=1\right)} \\=& \frac{\pi_{k} \mathcal{N}\left(\mathbf{x} | \boldsymbol{\mu}_{k}, \mathbf{\Sigma}_{k}\right)}{\sum_{j=1}^{K} \pi_{j} \mathcal{N}\left(\mathbf{x} | \boldsymbol{\mu}_{j}, \mathbf{\Sigma}_{j}\right)} \end{aligned}$$

9.2.1 Maximum likelihood

假设数据集$\left\{\mathbf{x}_{1}, \ldots, \mathbf{x}_{N}\right\}$由 $N×D$ 的矩阵$\mathbf{X}$，其中第$n$行表示为 $\mathbf{x}_{n}^{\top}$，相似地，隐变量由 $N×K$ 的矩阵$\mathbf{Z}$，其中第$n$行表示为 $\mathbf{z}_{n}^{\top}$，如果数据满足独立同分布，那么此时概率图表示为：

按照式子(5)描述的GMM，其对数似然函数可以表示为：

$$\ln p(\mathbf{X} | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})=\sum_{n=1}^{N} \ln \left\{\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\mathbf{x}_{n} | \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)\right\}$$

直接最大似然以上函数存在多个问题：

• 对于多元高斯，可能其中一个高斯坍缩到一个点，方差无穷小，导致似然函数最大
• 解出一组参数后，不知那个对应1号，那个对应2号，即不可确认性
• 对数没有直接作用于概率，而是作用于求和公式，此时求导得不到闭环解

由于以上问题的存在，上面似然只能用非常困难的，约束很多的基于梯度的优化方法，反正难整。

9.2.2 EM for Gaussian mixtures

对于带隐变量模型的最大似然求解，一个优雅而强大的方法就是EM(expectation-maximization)算法，EM非常强大，这里首先展示一下EM如何求解GMM的。

首先把似然函数对$\boldsymbol{\mu}_{k}$求导：

$$\begin{array}{l}{\frac{\partial}{\partial \boldsymbol{\mu}_{k}}\left(\sum_{n=1}^{N} \ln \left\{\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\mathbf{x}_{n} | \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)\right\}\right)} \\ {=\sum_{i=1}^{N}\left(\frac{\pi_{k} \mathcal{N}\left(\mathbf{x} | \boldsymbol{\mu}_{k}, \mathbf{\Sigma}_{k}\right)}{\sum_{j=1}^{K} \pi_{j} \mathcal{N}\left(\mathbf{x} | \boldsymbol{\mu}_{j}, \mathbf{\Sigma}_{j}\right)} \right) \frac{\partial}{\partial \boldsymbol{\mu}_{k}}\left(-\frac{\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right)^{\top} \boldsymbol{\Sigma}_{k}^{-1}\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right)}{2}\right)} \\ {=\sum_{i=1}^{N} \gamma\left(z_{n k}\right) \frac{\partial}{\partial \boldsymbol{\mu}_{k}}\left(-\frac{\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right)^{\top} \boldsymbol{\Sigma}_{k}^{-1}\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right)}{2}\right)} \\ {=\sum_{i=1}^{N} \gamma\left(z_{n k}\right) \boldsymbol{\Sigma}_{k}^{-1}\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right)} =0 \end{array}$$

使其导数等于0，其中协方差对于N是常数，因此可以约掉，得到：

$$\boldsymbol{\mu}_{k}=\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma\left(z_{n k}\right) \mathbf{x}_{n}$$

其中$N_{k}=\sum_{n=1}^{N} \gamma\left(z_{n k}\right)$，在$K$-means 中，可以理解为分配到 $k$ 簇的点的数量

同理，对协方差矩阵 $\boldsymbol{\Sigma}_{k}$ 求导并置零，可以解出闭环解：

$$\boldsymbol{\Sigma}_{k}=\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma\left(z_{n k}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right)^{\mathrm{T}}$$

注意以上情况并没有考虑混合系数加起来等于1的约束，因此最后再求解混合系数 $ \pi_{k}$ , 用拉格朗日算子把系数约束加入最大似然：

$$\ln p(\mathbf{X} | \boldsymbol{\pi}, \boldsymbol{\mu}, \mathbf{\Sigma})+\lambda\left(\sum_{k=1}^{K} \pi_{k}-1\right)$$

对 $ \pi_{k}$ 求导并置零，得到如下，可见前项又是后验概率：

$$0=\sum_{n=1}^{N} \frac{\mathcal{N}\left(\mathbf{x}_{n} | \boldsymbol{\mu}_{k}, \mathbf{\Sigma}_{k}\right)}{\sum_{j} \pi_{j} \mathcal{N}\left(\mathbf{x}_{n} | \boldsymbol{\mu}_{j}, \mathbf{\Sigma}_{j}\right)}+\lambda$$

上式根据参数和约束，可以解出$\lambda=-N$ ，且混合系数为：

$$\pi_{k}=\frac{N_{k}}{N}$$

可以看出，式子(12)(13)(17)依次更新了均值，协方差和混合系数，不过由于后验参数$\gamma (z_{n k})$ 的存在，上面三个系数只能和后验参数之间迭代优化，首先随机初始化均值、协方差和混合系数，然后开始迭代更新，这就是GMM的EM方法，用当前参数值求解后验$p(\mathbf{z}| \mathbf{x})$的过程，叫做E-step, 用后验概率带入似然，再最大似然求解参数，这个过程叫做M-step。

总结一下EM for Gaussian Mixtures:

1. 初始化均值 $\mu_{k}$, 协方差 $\Sigma_{k}$ 和混合系数 $\pi_{k}$, 并计算似然函数初始值
2. E step. 根据参数求解后验概率，即 $p(\mathbf{z}| \mathbf{x})$ ：

$${\qquad \gamma\left(z_{n k}\right)=\frac{\pi_{k} \mathcal{N}\left(\mathbf{x}_{n} | \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)}{\sum_{j=1}^{K} \pi_{j} \mathcal{N}\left(\mathbf{x}_{n} | \boldsymbol{\mu}_{j}, \mathbf{\Sigma}_{j}\right)}}$$

1. M step.根据后验概率更新模型参数：

$${\mu}_{k}^{\mathrm{new}} =\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma\left(z_{n k}\right) \mathbf{x}_{n}$$ $${\Sigma}_{k}^{\mathrm{new}} =\frac{1}{N_{k}} \sum_{n=1}^{N} \gamma\left(z_{n k}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}^{\mathrm{new}}\right)\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}^{\mathrm{new}}\right)^{\mathrm{T}}$$ $$\pi_{k}^{\mathrm{new}} =\frac{N_{k}}{N}$$

​ 其中

$$\qquad N_{k}=\sum_{n=1}^{N} \gamma\left(z_{n k}\right)$$

1. 计算似然函数值：

$$\ln p(\mathbf{X} | \boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Sigma})=\sum_{n=1}^{N} \ln \left\{\sum_{k=1}^{K} \pi_{k} \mathcal{N}\left(\mathbf{x}_{n} | \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k}\right)\right\}$$

​ 根据似然的收敛情况继续迭代Step2。

9.3 An alternative view of EM

EM 算法的目标是找到带隐变量模型的最大似然解，假设数据集$\left\{\mathbf{x}_{1}, \ldots, \mathbf{x}_{N}\right\}$由 $N×D$ 的矩阵$\mathbf{X}$，其中第$n$行表示为 $\mathbf{x}_{n}^{\top}$，相似地，隐变量由 $N×K$ 的矩阵$\mathbf{Z}$，其中第$n$行表示为 $\mathbf{z}_{n}^{\top}$，模型参数用 $\boldsymbol{\theta}$ 表示，那么对数似然函数可以表示为：

$$\mathcal{L}( \boldsymbol{\theta})=\ln p(\mathbf{X} | \boldsymbol{\theta})$$

由于似然函数是带有隐变量的，并不是一个数值，因此最大化 $\mathcal{L}$ 实际应该是最大化$\mathbb{E}_{\mathbf{z}}[\mathcal{L}]$ ，所以：

$$\mathbb{E}_{\mathbf{z}}[\mathcal{L}]=\mathbb{E}_{\mathbf{z}}[\ln p(\mathbf{X,Z} | \boldsymbol{\theta})]$$

上式中，带隐变量的似然也就是把联合分布对隐变量求期望。对于每个观测变量 $\mathbf{X}$ 都有一个对应的隐变量 $\mathbf{Z}$ ，也就是说 $\{\mathbf{X}, \mathbf{Z}\}$ 在一起是一组完整的数据，称作complete 数据，对于单独一个 $\mathbf{X}$ 是不完全数据，称作incomplete 数据。对于一般问题，我们掌握的往往是不完全数据 $\mathbf{X}$ 而不是完全数据 $\{\mathbf{X}, \mathbf{Z}\}$ ，而我们对隐变量的认知一般表示为后验概率$p(\mathbf{Z} |\mathbf{X}, \boldsymbol{\theta})$，因此我们往往先根据观测值和初始参数，求解隐变量的后验概率，后验得到之后，期望也就可以表示出来了，因此这一步就是E-step, 后面的M-step就是最大化期望求解新的 $\boldsymbol{\theta}$。

在E-step中，用当前的参数 $\boldsymbol{\theta}^{old}$ 来求解后验 $p(\mathbf{Z} |\mathbf{X}, \boldsymbol{\theta}^{old})$ ，从而表示完全数据的期望，有了隐变量的分布，就可以用积分求期望了，此处引入目标期望的新表示$\mathcal{Q}\left(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {old }}\right)$：

$$\begin{aligned} \mathcal{Q}\left(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {old }}\right) =&\mathbb{E}_{\mathbf{z}}[\mathcal{L}] = \int_{z} p(\mathbf{Z})\mathcal{L} d z \\=& \int_{z} p(\mathbf{Z})\ln p(\mathbf{X,Z} | \boldsymbol{\theta}) d z \\=& \int_{z} p(\mathbf{Z} |\mathbf{X}, \boldsymbol{\theta}^{old})\ln p(\mathbf{X,Z} | \boldsymbol{\theta}) d z \\离散情况=& \sum_{\mathbf{Z}}p(\mathbf{Z} |\mathbf{X}, \boldsymbol{\theta}^{old})\ln p(\mathbf{X,Z} | \boldsymbol{\theta}) \end{aligned}$$

请注意，这个期望函数是人为设计成这样的，对数$\ln$ 直接作用在联合分布上，是为了最大化方便求解，其收敛性得到证明了。

在M-step中，通过最大化期望值来更新参数得到 $\boldsymbol{\theta}^{new}$ :

$$\boldsymbol{\theta}^{\text {new }}=\underset{\boldsymbol{\theta}}{\arg \max } \mathcal{Q}\left(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {old }}\right)$$

以上两步就是EM的通用更新方式，理论上，每次更新 $\boldsymbol{\theta}^{new}$ 都会使得$\mathbb{E}_{\mathbf{z}}[\mathcal{L}]$ 更大，直到收敛到极值。如果给了参数的先验知识 $p(\boldsymbol{\theta})$ ，EM算法也可以用于求解最大后验概率MAP，E-step是一样的，只是M-step加上了一项：$\mathcal{Q}\left(\boldsymbol{\theta}, \boldsymbol{\theta}^{\text {old }}\right)+\ln p(\boldsymbol{\theta})$ ，有了先验之后，GMM的奇异值问题就不会出现了。此外，EM算法也可以用于处理缺失值情况。贝叶斯EM，求解是变分法，就可以避免坍塌情况，因此，某些类别可以一个样本都没有，所以贝叶斯EM的K数量可以减小，有些K的概率始终是先验的概率，没有新的类别划过去。在普通EM中，会出现奇异值情况。

Down